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1. Introduction 

The phase structure of strongly interacting matter, i.e. the QCD phase diagram, has not yet 
been determined from first principles. The reason is of course well-known: at nonzero baryon 
chemical potential the fermion determinant is complex and this prevents the use of the standard 
numerical methods commonly used in lattice QCD simulations. As a result there is no agreed 
QCD phase diagram in the plane spanned by the temperature T and baryon chemical potential }1b 
(or quark chemical potential \l q = jUs/3). Some colourful impressions appearing close to the top 
of a Google search (yielding around 19000 hits) are shown in fig. |l[ 
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Figure 1: Some impressions of the QCD phase diagram. 



What is understood is the transition along the temperature axis for vanishing and small chem- 
ical potential. This is highly relevant since the crossover from the confined hadronic phase to the 
deconfined quark-gluon plasma phase is currently being probed by heavy-ion collisions at RHIC 
and the LHC. The theoretical status of the phase diagram in this region was reviewed by Maria- 
Paola Lombardo at this conference [1]. At larger chemical potential the sign problem prohibits 
the use of well-established methods. In this contribution I will first discuss the close connection 
between the sign, the overlap and the Silver Blaze problems and how they appear in numerical 
simulations. Subsequently I will present results in cases where simulations are possible, either 
because the sign problem is absent or because it is much milder than in full QCD. In the third and 
final part of this review I will motivate why it makes sense to go into the complex plane and discuss 
the applicability of complex Langevin dynamics. I will not review the evidence in favour or against 
the QCD critical endpoint, and quarkyonic and colour-superconducting phases. Also, "standard 
methods" at small jl/T, such as reweighting, Taylor series, analytical continuation, histograms or 
the canonical ensemble, have been partly reviewed at previous gatherings [2-4]. 

2. Sign, overlap and Silver Blaze problems 

In the standard formulation of lattice QCD, the quarks are integrated out analytically, yielding 
the fermion determinant, while the functional integral over the gluonic variables in the partition 
function remains to be done, i.e. 




DUD\ffD\tfe^ SyMiu) ^ U '^ 



J DUe- SYM< - u) detM(U;n). (2.1) 
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Here M is the fermion matrix, appearing in the femiionic action as 

S P (U;n) = - J d 4 xxj/M(U;ii)Y- (2.2) 

The chemical potential is introduced in the usual way, as the fourth component of an imaginary 
vector potential [5]; the dependence of M on the gauge links U will no longer be indicated. If 
the combined Boltzmann weight, from the Yang-Mills action and from the fermion determinant, 
is real and nonnegative, importance sampling can be used to evaluate the remaining bosonic high- 
dimensional integral numerically. This is the usual situation in lattice QCD at zero and nonzero 
temperature. 

At nonzero baryon density, or at nonzero baryon chemical potential in the grand-canonical 
ensemble, this is however not the case. Due to the lack of y 5 -hermiticity, 

Y 5 M(p)Y 5 =M\-tx*)^M\tx), (2.3) 

the determinant is complex and satisfies 

[detM(ju)]*=detM(-/i*). (2.4) 

The problems associated with a complex or nonpositive weight are collectively referred to as the 
sign problem, even though a complex phase problem would be more accurate. It is important to 
realize that the complexity does not arise from the Grassmann nature of the quark fields. After 
all, at vanishing chemical potential, the determinant is real. It is therefore not a fermion sign 
problem per se. In fact, the sign problem also appears in bosonic theories at nonzero chemical 
potential, in which the complex action enjoys the same symmetry properties as the determinant, 
i.e. S*(/J.) = 5(— ;U*). In all these cases the complexity is absent for purely imaginary chemical 
potentials. 

The simplest approach would be to write 

detM = |detM|e ! ' e , (2.5) 

and ignore the phase factor at first instance. Simulations in the resulting phase-quenched (pq) 
theory are straightforward in principle, since they do not suffer from the sign problem. The phase 
factor can be incorporated later via reweighting, 

_ [DUe- SYM detMO _ f DU e~ SYM | det M\e w O _ (e w 0) m 
^ '~ fDUe- s ™ det M ~ f DU e~ s ™\ det M\e ie ~ (e w ) pq ' 

so that observables are recovered from ratios in the phase-quenched theory. 

However, when the full and the phase-quenched theory differ, this is bound to lead to an 
overlap problem: the generated phase-quenched ensemble has little in common with the desired 
ensemble in the original theory. This is easily understood mathematically: the average phase factor 
( e ' d )pq m tne phase-quenched theory is nothing else than the ratio of two partition functions, 

e _ fDUe-s™\d*M\e» _ Z _ VAf/T 
{e )pq ~ fDUe-s™\ detM \ ~Z n ~ e ' /_/ /pq ' ( j 
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where f = — {T /V) InZ is the free energy density in the full theory, V the three-dimensional volume 
and T the temperature (and similar for / pq in the phase-quenched theory). When Af ^ 0, it follows 
that the average phase factor will go to zero in an exponential fashion when the thermodynamic 
limit is taken and the temperature is reduced. In lattice simulations at fixed lattice spacing, this 
scaling is with the four- volume, V /T ~ £1 = N^N t . This behaviour of the phase factor, and the 
related worsening of the signal-to-noise ratio, is unavoidable and limits the system sizes that can 
be addressed with this approach. If this happens, the sign problem is severe. Importantly, Af 
depends on the theory or model under consideration and in some cases, such as in QCD in the 
strong-coupling limit [6] or in certain effective three-dimensional models (see below), Af is smaller 
than perhaps expected and the sign problem can be rather mild. Simulations are then possible for 
reasonably large - but still finite - system sizes. 

There is also a physical reason for the presence of the sign problem. This is especially clear 
at zero temperature, in the case that the transition from vacuum to a high-density state (the onset) 
occurs at a lower critical chemical potential in the phase-quenched theory than in the original 
theory, i.e. jjic q < Mc- The complex phase should then produce excessive cancellations to wash out 
the early onset in the phase-quenched theoiy. This elimination has to be perfect, since in the original 
theory at zero temperature thermodynamic quantities are independent of jJ. as long as /J. < pL c . This 
has been named the Silver Blaze problem [7, 8] : while every numerically generated configuration in 
the ensemble depends on the chemical potential, this chemical potential dependence should cancel 
exactly in thermodynamic observables, when evaluated in the entire ensemble. The mechanism 
behind this might be highly nontrivial. The region where /i,P q < jU < jU c is referred to as the Silver 
Blaze region. 

A convenient example to illustrate this is QCD with two flavours of equal mass, and hence 
fermion determinant [detMQi)] 2 (jU is the quark chemical potential). The determinant in the phase- 
quenched theory is different and reads | detM(jU) | 2 = detM(/i) detM(— ju), i.e. the phase-quenched 
theory corresponds to a theory with nonzero isospin chemical potential [9]. The onsets occur 
at different critical chemical potentials: in the isospin theory one finds condensation of pions at 
jj.c q = m n /2, while at finite quark number one expects onset to nuclear matter at \i c ~ m^/3 (— 
binding energy), where mjy denotes the nucleon mass. The Silver Blaze region is therefore the 
region with m K /2 < /x < m^/3. It has indeed been shown that in the Silver Blaze region the jU- 
dependent spectral density of the Dirac operator is complex and highly oscillatory, with a diverging 
amplitude exp(£2) and period I/O, to ensure jU -independence of the chiral condensate and baryon 
number (OSV mechanism) [10, 11]. Only when the integration is done carefully can the desired 
cancellations take place. This important insight indicates some of the difficulties one may encounter 
when approximate methods are used. 

The discussion above refers to the case that the sign problem has not been resolved. However, 
in some cases the sign problem can be eliminated by a reformulation in other degrees of freedom 
or by the use of complex Langevin dynamics. In those cases, the distinction between mild and 
severe sign problems is no longer necessary. Similarly, there is no longer a Silver Blaze problem. 
However, it is still useful to study the approach in particular in the Silver Blaze region, since this 
provides a stringent test on whether the sign problem has indeed been eliminated correctly. 

After setting the stage, I will continue with some selected results in theories without a sign 
problem and subsequently discuss theories with a milder sign problem, in particular effective three - 
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dimensional models. Finally, I will argue that one has to take the complexity at nonzero chemical 
potential seriously and embark on a journey into the complex plane. 



3. Theories without a sign problem 

The prime example of a nonabelian gauge theory without a sign problem is two-colour QCD or 
QC2D: QCD with gauge group SU(2) and quarks transforming in the fundamental representation. 
In this theory, there is an additional Pauli-Giirsey symmetry [12], 



KM{pL)K- { =M*(H), 



K = Cy 5 x 2 , 



(3.1) 



with C the charge-conjugation matrix. While 75-hemiiticity is still absent at nonzero chemical 



potential, see Eq. ([2.3|), the Pauli-Giirsey symmetiy nevertheless leads to a real determinant, even 
when jU / 0. For two flavours and in the presence of diquark sources, the determinant can then be 
written in a manifestly positive manner and lattice simulations can be carried out using conventional 
algorithms. 
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Figure 2: Tentative phase diagram in QC2D [15]. 



Lattice QC2D at nonzero chemical potential has been studied already for some time. Recent 
progress concerns the quarkonium spectrum at low temperature and large chemical potential [13], a 
detailed analysis of singular values of the Dirac operator, determined by D' "(ju)D(ju) \\f n = ^\j/(n) 
[14] (in contrast to the eigenvalues of the Dirac operator, the singular values remain real and non- 
negative), and the phase diagram at nonzero temperature and density [15]. The latter is shown in 
Fig. ||, with the axes labelled both in MeV and in lattice units. Besides the usual confined and quark- 
gluon plasma phases, one encounters a low-temperature phase with a nonzero diquark condensate 
(qq) as well as a quarkyonic phase. While the transition from vacuum to a diquark-condensed 
phase should not appear in QCD with three colours, in QC2D it is indeed the expected behaviour, 
since baiyons are bosons and the lightest baryon is degenerate with the pion. The evidence for the 
quarkyonic phase is discussed in Ref. [16]. 
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However, QC2D is not the only theory with a real and positive determinant at nonzero chemical 
potential. Another interesting example is given by using as gauge group the exceptional group 
G2 [17]. Unlike QC2D, this theory has bosonic (qq) and fermionic (qqq) baryons, as well as 
mesonic (qq), hybrid (qggg) and multi-quark bound states. Thanks to an extended Pauli-Giirsey 
symmetry, the determinant is real and positive, even for a single flavour at nonzero /x. 



Finite-Li/T dependence || Finite-n/T dependence || Finite-bi/T dependence 
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Figure 3: Gauge group G2: chemical potential dependence of the Polyakov loop (left), the normalized chiral 
condensate (middle) and the normalized baryon density (right) for three temperatures [17]. 

Some first results are presented in Fig. |J where the Polyakov loop, the chiral condensate and 
the baryon density are shown as a function of chemical potential for three temperatures. The onset 
is presumably due to a diquark condensate, as in QC2D, but there are indications for a second rise 
in density at larger chemical potential (not shown). To settle this will require precise calculations 
of the spectrum at zero temperature. The drop of the Polyakov loop at larger /I is most likely a 
saturation effect, see e.g. the behaviour of the density on the right, and is hence a lattice artefact. 
Nevertheless, these first results are quite encouraging and clearly deserve further effort. 

4. Theories with a milder sign problem 

From the perspective of lattice QCD, effective lower-dimensional field theories or spin models 
provide a useful entry into the study of the QCD phase diagram for a number of reasons. First of 
all, it has become clear that these models may have a milder sign problem and are therefore more 
amenable to standard lattice techniques, in combination with reweighting. Secondly, sometimes 
they can be reformulated using different (or dual) degrees of freedom, such as in world line or flux 
representations, in which the sign problem is manifestly absent. And finally, some of these models 
can be solved numerically with complex Langevin dynamics, providing a completely independent 
approach. One role these models therefore have is to provide a useful testbed for new methods. 

Ideally, the effective models are constructed in such a way that a detailed mapping between 
the couplings appearing in those and in full QCD is available: the insight obtained can then be 
translated back to QCD. To illustrate this, I start with the simplest three-dimensional SU(3) spin 
model as a guide. The action, appearing as exp(— S), consists of two parts and is written as 

S = S B + S F , (4.1) 

with 

S B = - J8 £ [P X P; + P*P y ] , S F = [e*P x + e-l 1 !*] . (4.2) 

(xy) 
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Here the degrees of freedom are effective Polyakov loops, P x = TrU x , P* = Tr£/J, where the U x s 
are SU(3) matrices, living on a three-dimensional lattice. The QCD gauge action has been re- 
placed by a nearest-neighbour interaction, while static (anti)quarks are represented by (conjugate) 
Polyakov loops, weighted with the chemical potential to introduce an imbalance. The action is 
indeed complex and 5*(jU) = S(— jU*). The effective parameters can be mapped to the parameters 
in four-dimensional lattice QCD: at leading order in a combined strong-coupling and hopping ex- 
pansion, one finds that j8 ~ (Aid/ 18)^, h = (2k) n \ with K the Wilson hopping parameter, and jU 
corresponds to /J./T in four dimensions. I come back to this below. 

This model is of course already quite old: it was studied in a mean-field approximation nearly 
30 years ago [18] and it was one of the earliest models studied with complex Langevin dynamics 
[19,20]. Nevertheless, it has seen a recent revival for a number of reasons: it can be reformulated as 
a flux model without a sign problem [21]; the applicability of complex Langevin has been studied 
in detail [22] ; a renewed mean-field analysis is available [23] ; and in fact it is part of a whole family 
of high-order strong-coupling models, which have been studied in a systematic fashion in recent 
years [24]. 

Let us start with the flux representation, in which the sign problem is absent [21]. In this 
approach the Boltzmann weight is expanded, as in a classical high-temperature expansion, to all 
orders and the partition function takes the form of infinite sums of integrals over powers of (conju- 
gate) Polyakov loops at each site x, i.e., 

r(n x ,n x ) = ( dU x (TvU x ) n * (TrU*)** . (4.3) 

JSU(3) 

The crucial next step is to perform these single-site SU(3) integrals, which is possible in this case. 
The result is a monomer-dimer system with constraints, since I(n x ,n x ) is only nonzero provided 
(n x — h x ) mod 3=0. The important finding is that all the nonzero weights that appear in this 
representation are real and positive, even when ji^O. The resulting model can then be solved with 
importance sampling or a worm algorithm [25]. 



Endpoint for h = 0.005 




1 2 3 4 5 

Figure 4: SU(3) spin model: expected phase structure (left) and the result from simulations in the flux 
representation (right) [25]. 

This theoiy is expected to have a phase diagram as in Fig. ^ (left). For small fermion coupling 
h, a quasi-confined or disordered phase (with small (P) , recall that the centre-symmetry is explic- 
itly broken by the fermionic contribution) and a deconfined or ordered phase (with large (P)) are 
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separated by a first-order phase boundary in the /3 - /I plane. This line ends in a critical endpoint 
after which only a crossover remains. The simulations [25] confirm this, see Fig. |] (right), where 
phase boundaries are shown for various values of h. Increasing h, which corresponds to reducing 
the quark mass, lowers the value of j8 c and weakens the transition, as does increasing }i. Note that 
this is indeed the expected behaviour in full QCD with heavy quarks, as will be discussed below. 



The SU(3) spin model ( |4.2[ ) is in fact part of a whole family of effective Polyakov loop models 
constructed recently [24,26,27]. Systematically integrating out the spatial links in lattice QCD us- 
ing a strong-coupling expansion while leaving the temporal links untouched yields effective three- 
dimensional actions with more and more Polyakov loop interaction terms. For pure SU(3) gauge 
theory, these take the form [26] 

S = - Ai £ [P X P* + P*P y ] - A 2 £ [P X P* + P*P y ] + (higher-order representations) + . . . , (4.4) 

where nearest and next-to-nearest neighbour terms are indicated explicitly. Within each term, cer- 
tain strong-coupling contributions can be resummed to all orders, using e.g. 

s = -A! £ [p x p;+p;p y ] +... = -£ in [i + A! (p x p;+p x *p y )] +.... (4.5) 

(xy) (xy) 

Finally, the dependence of the couplings A,- on /3 and N T is known to quite high order; for instance 
for N? = 4 one finds [26] 



Xi(u,N x = 4) = M 4 exp 



( 4 ^ fi 1035317 , n \ 

4(4m 4 + 12 M 5 -14m 6 + ...+ gi2Q M 1Q + ...j 



(4.6) 



X 2 (u,N? =4) = w 8 [I2w 2 + 26w 4 + 364m 6 + ...] , (4.7) 

where u = w(j3) = a/(j3) = j8/18 + . . . , as it appears in the character expansion. By including 
higher-order terms, the approach can be systematically improved and errors due to the truncation 
can be estimated. 

Quarks are taken into account via a hopping expansion; for heavy (static) quarks this yields 
the well-known contribution [28-31] 



^lndet 



2 / ,„ ,\ 2 



1 + he^ /T W x 1(1+ he-^ /T 



(4.8) 



where W x is the untraced Polyakov loop, h = h{K : N z ) and the remaining determinant is in colour 
space. At lowest order in h and ff, with h = (2k) Nt , an expansion in h yields the fermionic con- 
tribution given in Eq. (JO). An important distinction between the two expressions is that the 



action in Eq. (4.8) is still a determinant and hence incorporates Fermi-Dirac statistics and satu- 
ration on the lattice, while these properties are lost in Eq. (fO|). This can e.g. be seen by writing 
down the expression for the density (n) = (T /V)dlnZ/dfi, using the SU(3) identity det (1 +cU) = 
l+cTr£/ + c 2 Tr£/ f + c 3 . 

It is clear that the approach sketched here goes considerably beyond the simplest SU(3) spin 
model discussed first and one may ask the question whether quantitative results can be obtained. 
In order to answer this, it may for a moment be advantageous to forget about the strong-coupling 
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Figure 5: Validity of strong-coupling models, as complementary to dimensional reduction (left). The critical 
temperature in SU(3) gauge theory from a one-coupling effective theory (right) [24]. 



origin and view the setup instead as a generic approach to construct effective models, whose valid- 
ity may be justified afterwards. In some sense it can be viewed as complementary to dimensional 
reduction at high temperature and weak coupling. This is illustrated in Fig. || (left). While dimen- 
sional reduction, formulated in terms of quarks and gluons, is formally only valid at weak coupling 
and hence (very) high temperature, it is well-known that it can be applicable down to \.5-2T c [32]. 
Similarly, while a combined strong-coupling/hopping expansion, formulated in terms of colour sin- 
glets and meson/baryon degrees of freedom, is formally valid in what appears to be an unphysical 
regime at finite lattice spacing (i.e. not in the continuum limit), one may still test its validity. This 
is illustrated in Fig. || (right) in pure SU(3) gauge theory. Here the critical temperature from four- 
dimensional lattice simulations (270 MeV) is compared with the one obtained from an effective 
three-dimensional theory with one coupling only, X\{fi,N x )- Once the critical coupling Ai iC (/$,iV T ) 
is determined, the results can be mapped back to four dimensions to find j3 c (Af T ). After setting 
the scale, this yields Fig. || (right), where all the data points come from a single simulation in the 
one-parameter effective theory. By comparing truncation effects, an error can be assigned as well. 

With the inclusion of heavy quarks, the top-right corner of the Columbia plot can also be 
studied. The Columbia plot, see Fig. || (left), indicates the expected quark mass dependence of the 
thermal deconfinement transition at vanishing chemical potential: a first order transition for light 
and heavy quarks, related to chiral and centre symmetry, and a crossover for intermediate-mass 
quarks. An extension to nonzero chemical potential makes the Columbia plot three-dimensional, 
as shown in Fig. [7] for the heavy quark corner: the vertical axis indicates (jx/T) 2 , with negative 
values corresponding to imaginary }X = ijX\. Due to the interplay with the centre-symmetry, a 
rich phase structure emerges, with a periodicity in the imaginary /I direction with period 2nT /3 
[33]. In the high-temperature deconfined phase this results in a Roberge-Weiss (RW) transition 
at jUrw = inT/3, where the expectation value of the Polyakov loop is rotated by an element of 
the centre and the trivial vacuum at vanishing /I is no longer preferred. The first-order transition 
at /J. = /Irw is only present at high temperature; in the confined phase at low temperature, the 
periodicity is smoothly realized. The RW transition at fixed pL\ therefore necessarily ends at an 
endpoint at a temperature Trw, see Fig. || (right) [33-35]. 
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Figure 6: Columbia plot: expected quark mass dependence of the order of the thermal deconfinement 
transition at vanishing chemical potential (left) and the extension of the QCD phase diagram to nonzero 
chemical potential in the case of heavy quarks (right). 




Figure 7: The extension of the top-right corner of the Columbia plot to nonzero chemical potential [24]. 



For light and heavy quarks, the thermal deconfinement transition is first order. It is expected 
that this transition is connected with the RW transition at imaginary ju, leading to the phase dia- 
gram in Fig. |6] (right). The RW endpoint is then a triple point, where three first-order lines meet. 
For quarks with an intermediate mass, the deconfinement transition is a crossover and the RW 
endpoint is a second-order point at the end of a first-order line. This has been confirmed with lat- 
tice QCD simulations in the two-flavour [36] and the three-flavour [37] case (the triple point has 
also been seen with holographic methods [38]). Returning to the Columbia plot in Fig. 0, these 
results imply that the plane at fixed (/J./T) 2 = — (ti/3) 2 is critical, with a tricritical line separating 
the first-order region for heavy quarks from the second-order region for intermediate-mass quarks. 
Moving away from the RW plane, the tricritical line turns into a critical (second-order) surface, 
which separates the thermal first-order and crossover transitions. An intriguing recent insight is 
that the tricritical line plays an important role in determining the curvature of this critical surface at 
(jlt/r) 2 > — (7l/3) 2 , via tricritical scaling [37]. The shrinking of the first-order region as (jx/T) 1 
is increased, according to tricritical scaling, has been confirmed in the effective Polyakov loop 
models discussed above [24] as well as in the three-state Potts model [37], simulated some time 
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ago [39]. Clearly, this goes far beyond simple analytical continuation and motivates the study of 
the three-dimensional Columbia plot using imaginary chemical potential also for light quarks [40]. 

5. Into the complex plane! 

For the remainder of this overview, let us take a few steps back and tiy to tackle the problem 
of a complex Boltzmann weight head on. Consider first the simple integral 

/oo 
dxe~ s{x) , S(x) = -ax 2 + ibx. (5. 1) 

-oo 2 

In Kindergarten one used to be taught to do this integral by completing the square (or by a saddle 
point approximation, for more complicated functions), with the result that the variable is taken into 
the complex plane, x — > x — ib/a. The lesson to be drawn from this is that one should not insist 
on staying on the real axis, but be more imaginative. This leads to a radically different approach: 
all degrees of freedom are complexified, x — >• z = x + iy, yielding an enlarged configuration space 
with literally new directions to explore. In SU(A0 gauge theories or matrix models, this extension 
implies that the dynamics takes place in SL(iV,C) [31,41]. The hope is that a real and positive 
distribution exists in this complexified space, as illustrated in Fig. ||. 




Figure 8: Weight in the partition function: while the original weight p(x) on the real axis is complex and 
oscillatory (left), there may be a real and positive weight P(x,y) in the complex plane (right). 

The question is of course how to find this distribution. Before discussing complex Langevin 
dynamics, I briefly mention a recent proposal for studying high-density QCD on a Lefschetz thim- 
ble [42]. The idea here is to deform the integration contour into the complex plane and associate 
(real) integration domains with each stationary point, the so-called thimbles. The original path 
integral is then expressed as a sum over integrals over thimbles. This proposal, partly motivated by 
ideas due to Witten and Morse theory, is currently being implemented and results, first for cases 
simpler than QCD, will be available soon [43]. 

One potential way to obtain a real and positive probability distribution P(x,y) in the complex- 
ified space is via the solution of a stochastic process. This is exactly the aim of complex Langevin 
(CL) dynamics, proposed many years ago by Parisi [44] and by Klauder [45] and the topic of the 
second half of this contribution. Since CL has been around for 30 years, the natural question to ask 
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is why to talk about it here. Moreover, while for real actions the procedure is equivalent to path 
integral quantization [46] and known as stochastic quantization [47], for complex actions there is 
no formal proof of correctness and "disasters of various degrees" were encountered by e.g. Amb- 
j0rn and others in the mid-1980s [48,49]. However, in recent years examples in field theory have 
been given for which CL can solve severe sign and Silver Blaze problems and (importantly!) also 
gives the correct result. In addition, our analytical understanding has improved steadily. In the 
literature various scattered results can be found, starting from the mid-1980s; here I will review 
mostly results obtained from 2008 onwards in the context of finite density [31] . 

5.1 Sign and Silver Blaze problems 

Let us first consider two examples where CL has successfully solved a severe sign problem 
as well as a Silver Blaze problem: the four-dimensional weakly interacting relativistic Bose gas 
[50,51] and one-dimensional QCD [52], both at nonzero chemical potential. 




Figure 9: Density («} in the relativistic Bose gas as a function of jj. in the phase-quenched (left) and the full 
theory (right), for four different lattice volumes (with am = X = 1). The onset occurs at \l c ~ 1 . 15 [50]. 

As in the case of QCD, the continuum and lattice actions for the Bose gas satisfy the usual 
symmetry, 5*(ju) = S(—fl*). For instance, the continuum action is 

S = J d A x[\d v ^\ 2 + {m 2 -^ 2 )\^\ 2 + ^{^*d^-d^*^) + X\^\ A ]. (5.2) 

At zero temperature, the theory has a Silver Blaze problem: in the full theory onset occurs at 
jj, c ~ m, with jU -independence when jx < /i c . On the other hand, in the phase-quenched theory the 
term linear in /I is absent, leaving just an effective /I -dependent mass term, m 2 s = m — \l 2 , and 
there is immediate \i -dependence as \i is increased. The purely imaginary term, linear in \i, should 
cancel this \i -dependence in the full theory (the same considerations hold of course for the lattice 
discretized theory) [50]. This is illustrated in Fig. ^. Shown are the densities (n) in the phase- 
quenched (left) and the full (right) theory, obtained with real resp. complex Langevin dynamics, on 
four lattice volumes from 4 4 to 10 4 (lattice units are used, with m = X = 1). In the full theory, the 
density is approaching zero in the thermodynamic and zero-temperature limit, until the transition 
to a condensed phase occurs at /I ~ 1.15, in agreement with mean-field expectations [51]. In the 
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Figure 10: Average phase factor in the phase-quenched relativistic Bose gas for four different lattice vol- 
umes, indicating the severity of the sign problem [50]. The lines are mean-field predictions (left). Difference 
in free energy densities between the full and the phase-quenched theory, in the mean-field approximation [5 1 ] 
(right). 



phase-quenched theory the density increases immediately, again in agreement with mean-field ex- 
pectations. The effects of the complex phase, e' 6 = e~ s /|e~ s |, are therefore correctly incorporated 
by complex Langevin dynamics. Furthermore the thermodynamic limit poses no problem. 

The severity of the sign problem can be studied by analyzing the average phase factor in the 
phase-quenched theory. Recall that 

(j e )n = ^- = e-°* f , A/ = /-/„, Q = #X- (5-3) 

The result is shown in Fig. [l0| (left). It can be seen that the phase factor behaves exactly as expected: 
it is 1 when jJ. = 0, but as soon as /i > it decreases, in agreement with the density plots above. 
As the volume is increased at fixed fi, the phase factor drops exponentially. Around the critical 
chemical potential, fi c ~ 1.15, the phase factor vanishes even on the 4 4 volume; yet this does not 
affect the detection of the phase transition. It is possible to compute the difference in free energies 
Af in mean-field theory on the lattice [51]. The result is shown in Fig. [T^ (right) in the symmetric 



phase. Using this in Eq. (5.3) yields the curves in Fig. 10 (left), in agreement with the numerical 
data. We conclude that CL can solve this theory, with a severe sign and Silver Blaze problem, in 
the thermodynamic limit. 

One-dimensional QCD [53-55] is exactly solvable and shares many features with QCD in the 
Silver Blaze region: the phase-quenched theory has a transition at pL = pL c = arcsinm, with m the 
quark mass, while the full theory has no transition at all, at least when the gauge group is U(iV c ). 
The sign problem is severe when jj. > ;U C . The chiral condensate can be written as an integral over 
the spectral density, 

L=[d 2 z^l, (5.4) 
J z + m 

where p(z;n) is complex, oscillatory and /I -dependent [55]. The condensate itself is /i -independent 
(Silver Blaze): for this to happen, the integral ( |5.4[ ) has to be carried out with great precision. When 
/I > jU c , this is nontrivial: after writing z = (e' x+ ^ — e _w ~^)/2, the spectral density has a diverging 
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Figure 11: One-dimensional QCD: Chiral condensate as a function of the quark mass for two lattice vol- 
umes n (left). Complex Langevin solution to the sign and Silver Blaze problems: the original complex 
and oscillatory weight p(x) is replaced by a fiat and localized distribution P(x,y) in the complex plane 
(right) [52]. 

amplitude and vanishing period in the thermodynamic limit, p ~ e n(n-Hc)+mx^ w here n is the one- 
dimensional volume. Only when all the oscillations are integrated precisely does the /i -dependence 
disappear. CL has no problems with this: the result for the chiral condensate is shown in Fig. [II] 
(left) as a function of ji c at fixed }i = 1 and is seen to agree with the exact result. The discontinuity 
at jU c = m = emerges in the thermodynamic limit. 

In this case, the distribution P(x,y) in the complexified space can be computed analytically 
in the thermodynamic limit. It is surprisingly simple: a flat distribution which is nonzero only 
in a strip in the complexified space, ji. — ji. c < y < ji. + ji c , see Fig. [II] (right). The complicated 
dynamics on the real axis has been replaced by free diffusion in the complex plane and CL finds 
this distribution as a solution of the stochastic process, exactly as it says on the tin [52]. 

5.2 Analytical insight 

As mentioned above, complex Langevin dynamics has a troubled past. Problems can be di- 
vided into three groups [56]: (i) there can be numerical runaways and instabilities, (ii) the theo- 
retical status is unclear, (iii) there may be convergence to the wrong limit. The first problem is 
numerical and follows from properties of the classical drift terms appearing in the CL equation. 
It can be partly eliminated by choosing a small enough Langevin stepsize [57] and solved by in- 
troducing an adaptive stepsize [58]. This has been implemented in a variety of theories, including 
SU(3) gauge theory, and no instabilities have been observed. It should be mentioned that theories 
which suffer from extreme instabilities and runaways are also prone to converge to the wrong limit 
(see below), so it is likely that the two are related. 

The statement that the theoretical status is unclear originates from the following: for real 
actions it is relatively straightforward to show that stochastic quantization is equivalent to path 
integral quantization [46]. This relies on the appearance of a semi -positive Fokker-Planck hamil- 
tonian. Unfortunately, this path can no longer be followed when the action is complex, since the 
associated Fokker-Planck hamiltonian is now no longer semi-positive. Recently, the theoretical 
foundation has been clarified and necessary conditions for obtaining the correct result have been 
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identified [56, 59]. The idea is as follows: let us consider again one degree of freedom x and the 
weight p(x) ~ e~ s M with S(x) complex. What is wanted are expectation values of the form 

(0{x,t)) p = J dxp(x,t)0(x), (5.5) 

where p(x,t) satisfies a (complex- valued) Fokker-Planck equation (t is the Langevin time) 

d r p(x,t) = d x (d x + S'(x))p(x,t), (5.6) 

with stationary solution p (x) ~ . Instead, what is solved in the CL process are expectation 
values of the form 

(0(x,t)) P = J dxdyP(x,y,t)0(x + iy), (5.7) 
where P(x,y,t) satisfies a (real-valued) Fokker-Planck equation 

d t P(x,y,t) = [d x (dx - K x ) - d y K y ] P(x,y,t), (5.8) 

with^ r = — ReS', K y = — ImS'. The crucial question is whether (0(x,t))p = (0(x,t)) p , possibly for 
all times or perhaps in the limit t — > oo. The answer, which relies on the use of the Cauchy-Riemann 
equations and partial integration, turns out to be yes, provided some conditions are met [56,59]: 
the important one is that the distribution P(x,y) drops off rapidly in the imaginary y direction, such 
that partial integration without boundaiy conditions is possible. A complication here is that the FP 
equation ( |5.8| ) cannot easily be solved, even in simple models (some recent examples can be found 
in Refs. [51, 56, 60]). A way forward is to require that partial integration is possible for products 
of the form 0(x + iy)P(x,y), for a large enough set of operators 0(x). This condition can then be 
translated into criteria for correctness, which can be written as (LO(x + iy)) = 0, with L a Langevin 
operator, and measured during a numerical simulation (see below). In any case, the main lesson is 
that the distribution should be sufficiently localized. In the case of one-dimensional QCD discussed 
above, this is achieved perfectly since P(x,y) vanishes outside a strip. 

The third problem of convergence to the wrong limit is the hardest one. It had been observed 
in the early days [48,49] and a more recent study in the three-dimensional XY model at nonzero 
chemical potential revealed an intriguing connection with the phase structure in the jS — jU plane 
[61]. I will come back to this and possible resolutions towards the end. 

5.3 The SU(3) spin model revisited 



I now return to the SU(3) spin model of Eq. ( |4.2| ) to illustrate some of issues discussed above. 



While this model had been solved with CL earlier [19,20], no detailed justification was provided. 
This was rectified in Ref. [22], in which the reliability of CL was studied using a variety of tests. 
Here I will present some results based on analyticity in p 2 , a low-order Taylor expansion, the 
criteria for correctness mentioned above, and finally a comparison with the flux representation. 

In Fig. H (left) (Tr (U + U^)/2) is shown as a function of p 2 (note that Tr{7 and Tr^" 1 
were earlier denoted with P and P*). At nonzero real p, (TrU) and (TrU~ l ) are both real but not 
identical, while at nonzero imaginary p they are complex with opposite imaginaiy parts. The sum 
is always real and, since it is even under charge conjugation, a function of p 2 instead of p. For 
imaginary p the action is real and the data has been obtained with real Langevin dynamics, for real 
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Figure 12: Analyticity in /I 2 in the SU(3) spin model: (Tr(I/ + £/~ 1 )/2) as a function of /I 2 for a number 
of j3 values. Data at imaginary fi (with fi 2 < 0) has been obtained with real Langevin dynamics, data at real 
jU (with fi 2 > 0) with complex Langevin dynamics (left). (Conjugate) Polyakov loops (TrU) and (TrU^ 1 ) 
as a function of ji. The inset shows a comparison with the lowest-order Taylor expansion (right) [22]. 



jj. complex Langevin dynamics has been used. The jS values are chosen both below and above the 
critical j6 f ~ 0.133 at /i = 0. The transition gets weaker as /I 2 is increased, in agreement with the 



Columbia plot discussion for heavy quarks. The important message from Fig. 12 (left) is that the 
numerical data at jj. > appear to be completely consistent with those at negative jj. 2 , where the 
reliability of Langevin dynamics poses no problem. This is a first test CL should pass (and which 
it fails in e.g. the symmetric phase of the three-dimensional XY model [61]). 

As a second test we study the splitting between (TrU) and (Trf7~ ! ) at small /». This model 
does not have a Silver Blaze region, as can be seen from e.g. the expression for the density 
(n) = (he^TrU — he^^TrU^ 1 ). It follows from the absence of the possibility to take the zero- 
temperature limit, since N x is no longer present as a free parameter. At nonzero pL, (TrU) and 
(TrU~ l ) differ therefore immediately and behave as 

(Tr?7) = c 1 +c 2 /2 J u + ^( J u 2 ), (TrU- 1 ) = c l -c 2 hll + 0(lJ, 2 ), (5.9) 

where the coefficient of the linear term c 2 is determined from a Polyakov loop correlator (note that 

c 2 <0) 

This is illustrated in Fig. |l2| (right). The inset shows that the CL results agree with the lines 
predicted by the first-order Taylor expansion at very small }i. 

The tests described so far use the connection with vanishing or imaginary /I and are therefore 
only applicable for small values of fi. In order to assess larger /J. values, we turn to the criteria for 
correctness, (LO) = 0, with L a Langevin operator and O = TtU,TtU~ . The results are presented 
in Fig. 13: the figure on the left shows (TrU) and (TrU^ 1 ) as a function of the Langevin stepsize 
e, while in the figure on the right the criteria are shown. Some explanation is in order. Numerical 
solutions of discretized versions of the Langevin equation, with stepsize £, suffer from finite step- 
size corrections. For the lowest-order standard discretization, these corrections are linear in e. This 



is clearly visible in Fig. 13, both for the observables and the criteria (red symbols). It is possible 
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Figure 13: Langevin stepsize dependence of (Trt/) and (TrU ! ) (left) and the real part of (LTrt/) and 
(LTrC/ -1 ) (right) using the standard lowest-order (red) and the improved (blue) algorithm [22]. 



to improve this behaviour by implementing a higher-order algorithm. In Ref. [22] we have chosen 
for the algorithm described in Ref. [62] (for real Langevin dynamics), which is easy to implement. 
Analytically, one expects a cancellation of the linear term in e, but a remaining dependence on £ 3 / 2 , 
but we found effective stepsize independence in the case of the observables - see the blue symbols 



in Fig. |13| (left) - and a much reduced stepsize dependence for the criteria - blue symbols in Fig. |13| 
(right). In the limit that e — > 0, both calculations agree and the criteria for correctness go to zero, 
as it should be. As a general remark, we found that the criteria for correctness are very sensitive: 
for instance also in the case of real actions and real Langevin dynamics, they indicate the presence 



of finite stepsize corrections [63]. We also note that the chemical potential in Fig. 13 is not small 
(jit = 3): the analyticity around jU = is then no longer of use and a justification as discussed here 
is necessary. 
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" 0.130 
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Figure 14: (Tr(I7 + U~ l )/2) as a function of jx 2 , as obtained with the flux algorithm (filled symbols), the 
lowest-order CL algorithm (empty symbols) and the improved CL algorithm (crosses). The values at jj. = 
(asterisks) have been obtained with importance sampling in the original formulation [25]. 



A final justification comes from a comparison with the flux representation discussed above. 
In Fig. 14 the \i 2 > side of Fig. 12 (left) is shown again, but in this case there are also symbols 
representing the results obtained with the flux formalism, which is sign-problem free. Overall, 
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agreement is seen, except in the critical region at j8 = 0.132. Interestingly, the discrepancy is also 
present at /i = 0, where the action is real, which suggests that it is not an issue related specifically to 
complex Langevin dynamics. Instead, the problem lies with the finite Langevin stepsize, since the 
Langevin results in Figs. |l2| and 14 have been obtained using the lowest-order algorithm without 
an extrapolation to zero stepsize. This can be remedied with the improved algorithm, and indeed 
the data indicated with the blue crosses is seen to agree with flux representation. This agreement is 
also present at larger jU = 3 and in both the ordered and disordered phases [25]. The conclusion is 
that CL is seen to work in the three-dimensional SU(3) spin model in the entire phase diagram. 



5.4 Modifications of the complex Langevin process 

In the SU(3) spin model, CL was seen to pass all the tests. The important question is to 
understand why. The crucial observation here is that the distribution in the complexified field 
space turns out to be sufficiently localized and that there is fast decay in the imaginary direction, as 



required by the mathematical justification discussed in Sec. p.2| . In contrast, in the disordered phase 
of the XY model the distribution is broad - there is slow decay - and CL fails [61]. It turns out 
that the relevant feature that distinguishes the XY model from the SU(3) model is that the latter is 
nonabelian, with a nontrivial Haar measure. It is exactly the presence of the reduced Haar measure 
in the spin model that ensures that the complexified field space is explored in a controlled fashion, 
with a localized distribution, and that the real manifold is stable under small fluctuations, since its 
contribution to the drift is purely restoring [64]. 

The question is whether this insight can be employed in a practical manner to improve wrongly 
converging processes. Since the reduced Haar measure is essentially a Jacobian, one may con- 
sider generating a nontrivial Jacobian by a field redefinition, chosen in such a way that the re- 
sulting force is restoring in the complexified space. Consider an action S{x) with a drift term 
K{x) = —S'(x), which happens to lead to unstable or wrongly converging dynamics after com- 
plexification, x — > x + iy. A field redefinition, x = x(u), will then lead to a new process, with the 
action S e s(u) = S(x(u)) — \nJ{u), where J(u) = dx{u)/du is the Jacobian. The corresponding drift, 
K(u) = —S' eii (u) = —S'(u) +J'(u) /J(u), is singular at J{u) = 0, but, importantly, may lead to a very 
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Figure 16: Complex Langevin results for the correlator (x 2 ) in the Gaussian model with a = a + i as 
a function of a after the transformation x = m 3 . The lines indicate the expected, analytically continued 
results [64]. 



different Langevin process after complexification, u — > u + iv. 

A simple and amusing example to illustrate how this works is the Gaussian integral, 

dxe-^ ax \ o = a + ib, (x 2 ) = — = ""'f, , (5.11) 

a a 1 + b l 

but with a negative real coefficient a < 0. After complexification, the classical drift is highly 



unstable as illustrated in Fig. |15| (left) and the complex Langevin process does not converge. A 
simple redefinition x(u) = u 3 changes the flow diagram dramatically. After complexification u — > 
u + iv, the contribution to the drift from the Jacobian, J'(u)/J(u) = 2/u, is always directed towards 
the real axis, stabilizing the flow, and new attractive and repulsive fixed points appear, as shown in 



Fig. [B] (right) [64]. 

Perhaps surprisingly, a CL computation of (x 2 ) = (u 6 ) in the new variables yields the cor- 



rect result (in the sense of being the analytically continued answer to negative a), see Fig. |16[ A 
field redefinition yields therefore a process which is manifestly distinct from the original one. As 
discussed in Ref. [64], this setup turns out to be related to the introduction of kernels in CL [65]. 

It is possible to take these ideas a step further in S\J(N) gauge theories and modify the Langevin 
process by interspersing dynamical updates with SL(iV,C) gauge cooling [66], with the aim of 
controlling the distribution in the extended field space by keeping the links as close as possible 
to being unitary. Exploiting the gauge freedom was earlier used in a slightly different way in an 
SU(2) toy model in Ref. [67]. Very recent results in QCD with heavy quarks indicate that this is 
a successful strategy, allowing for the first time access to high densities in SU(3) gauge theory, 
although there are still some open questions about the applicability at smaller j8 values [66]. 

6. Summary and outlook 

From what is described above, I conclude optimistically that the heavy-quark corner will be 
settled, using a combination of techniques, formulations and effective models. As we have seen in 
the past years, the interaction between conventional techniques, approaches based on new formula- 
tions which are free of the sign problem, and complex Langevin dynamics has been very stimulating 
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and resulted in new insights and solutions to old problems. On the other hand, in my view light 
quarks at nonzero density remain a challenge. The reason is that there is a strong sensitivity to 
quark masses and that effective models are much less predictive. 

For the case of complex Langevin dynamics, we have seen that the method can handle severe 
sign and Silver Blaze problems, phase transitions and the thermodynamic limit. However, a correct 
outcome is not guaranteed and there may be convergence to the wrong result. In the past years we 
have developed a better understanding, both on a formal level and with regards to numerical simu- 
lations, and developed a series of tests to justify (or dismiss) the results from numerical solutions. 
Importantly, this theoretical framework is also applicable to SU(3) lattice theory. 

Very recently, several proposals to modify the complex Langevin process have been put for- 
ward, of which gauge cooling [66] is the most promising. It will be interesting to investigate this 
further. 
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